Density Matrix Renormalization Group in the Heisenberg Picture 
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In some cases the state of a quantum system with a large number of subsystems can be approxi- 
mated efflciently by the density matrix renormalization group, which makes use of redundancies in 
the description of the state. Here we show that the achievable efflciency can be much better when 
performing density matrix renormalization group calculations in the Heisenberg picture, as only the 
observable of interest but not the entire state is considered. In some non-trivial cases, this approach 
can even be exact for finite bond dimensions. 
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Introduction - Quantum many-particle systems give 
rise to a number of intriguing phenomena such as 
quantum phase transitions, magnetic frustration, the 
existence of rare-earth magnetic insulators or high- 
temperature superconductivity. But as the size of the 
Hilbert space grows exponentially with the number of 
subsystems, the numerical simulation of such quantum 
many-body systems is difficult and often intractable. 

In some cases, however, a quantum system does not ex- 
plore its entire Hilbert space and numerical approaches 
like the density-matrix renormalization group (DMRG) 
technique [ij, become efficient tools. DMRG can be un- 
derstood as a variation over the set of matrix product 
states (MPS) @, S] whose size grows only polynomially 
with the number of subsystems. Its success is linked to 
the existence of an upper bound for the entanglement of 
contiguous sub-blocks of the system under study 0, 0, Q ■ 
This approach is therefore expected to work particularly 
well for the ground state of one-dimensional gapped sys- 
tems, in which correlation functions decay exponentially 
and the entanglement entropy saturates, satisfying an 
"area law" Q . There are of course situations in which 
no upper bound to the entanglement in the system exists 
or where it grows in time. In such cases the performance 
of DMRG deteriorates. This is typically the case for the 
dynamics of non-equilibrium states, as exemplified in re- 
cent studies of sudden quenches to Bose- Hubbard Hamil- 
tonians Due to its dynamical production, the entan- 
glement per unit area may grow linearly in time in those 
scenarios ■ To achieve a fixed precision DMRG al- 

gorithms hence need to use matrix dimensions that grow 
exponentially in time rendering them inefficient [l^ . 

An increasing number of experimental settings, includ- 
ing arrays of Josephson junctions ultra cold atoms in 
optical lattices . ion traps |l3| and arrays of coupled 
micro-cavities , offer the possibility to generate effec- 
tive many-particle systems. Hence, dynamical studies of 
quantum many-particle systems are expected to receive 
increasing attention in the future. Moreover in real ex- 



perimental situations, such systems will typically suffer 
from decoherence and dissipation and hence evolve into 
mixed states whose numerical description is even more 
demanding. It is therefore desirable to develop new more 
efficient methods for such problems or alternatively to 
improve existing DMRG methods further. 

In this letter we describe an approach to enhance the 
performance of DMRG in time-dependent settings by 
computing directly the observable of interest and prop- 
agating its evolution in the Heisenberg picture. Our ap- 
proach avoids calculating components of the quantum 
states that are irrelevant to the observable of interest and 
hence contrasts the standard DMRG approach, which 
generates a quantum state for the entire quantum sys- 
tem in the Schrodinger picture. DMRG performed in 
the Heisenberg and Schrodinger picture are thus comple- 
mentary. Whereas Schrodinger picture simulations cal- 
culate the entire state and subsequently allow to com- 
pute any observable of interest, Heisenberg picture cal- 
culations only consider one operator whose expectation 
value can in turn be computed for any initial state. 

In the sequel, we demonstrate that DMRG performed 
in the Heisenberg picture (H-DMRG) can have signifi- 
cant advantages for numerical simulations of quantum 
many-particle dynamics. These advantages become most 
significant in open system dynamics described by mixed 
states but can also be demonstrated rigorously for certain 
exactly solvable systems. We find numerical indications 
for a saturation of the block entanglement in the Heisen- 
berg picture for increasing system size which suggest that 
H-DMRG has superior efficiency in many cases. 

Main part - For linear chains of interacting sub- 
systems, we consider the evolution of operators such as 
X^{t) = U{t) (1®™-! ®X® 1®^-™) U{t)\ where X is 
a Hermitian operator acting on site m, and use a matrix- 
product representation 



Pi, 
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with suitable d x d-dimensional matrices A,-'^ and the 
canonical operator basis {Pq, Pi, P2, P3} with {Pm)ij — 
Sm,2i+j for i,j G {0, 1}. We focus our study on dynamics 
of the anisotropic Hcisenbcrg Hamiltonian for a chain of 
N spins, 



H 



N 



(2) 



as this model is known to exhibit dynamics that is nu- 
merically hard to simulate. In eq. Bz is an applied 
magnetic field, Jx, Jy and are spin-spin couplings and 
(tJ, (tJ and ctJ the Pauli operators at site j. 

Exact results - It is noteworthy that the time evo- 
lution of certain operators can actually be represented 
exactly by a matrix product operator with fixed finite 
bond dimension. For Hamiltonians of the form of eq. Q 
with Jz = 0, all local operators that transform under the 
Jordan- Wigner transformation jisj into local fermionic 
operators remain exact matrix product operators with 
fixed finite dimension for all times. Examples of such op- 
erators arc CT^ whose time evolution is an exact matrix 
product operator for matrix dimension d = A and gener- 
ally any product of Pauli-operators with an even number 
of (7^ or (7^ operators and any number of cr^ operators. 

To see this, let us first define the fermionic annihila- 
tion and creation operators [l5l |. c„ 
icr-'„)/2. In terms of c,„ and 
(HD with = reads H = -S2^^^(2c]cj - I) + 

Jx Y!^=i{c]-Cj){c]^^+Cj+i) - Jy EpL^(c] +Cj)(4+i - 
Cj+i). Given that the Hamiltonian is quadratic in Cm 
and cj„, the Hcisenbcrg time evolution of an individual 
Hcisenbcrg operator such as Cm{t) is found to be Cm{t) ~ 

'Yli^^i + /3j(i)c]y In the fermionic picture this 

may be written as matrix product operator with matrices 
of dimension 2 as it is essentially the same as a W-state. 
Rewriting the rhs in terms of Pauli operators we find 

c^{t) = Ef=i («. it) nCi <^!4 + Pj (t) nCi ^i^j) , 

where = ^ (cr^' ± icr^). This in turn may be written 
as a matrix product operator of the form eq. ([1]) whose 



the Hamiltonian 



matrices have the structure A 



P-^^A. 



(N) 



aN{t)P2, 4'^ = Mt)Po, A 
i(i) 



(")-a™(t)P2,<' 



(m) 
2 



Pmit)P2, A 

(3^{t)P2, = -A, 4"'^ = ^0 - P3, ) = Po. 

As every spin operator may be expressed as a sum 
of products of fermionic operators we can now under- 
stand the above observations. For example, because 
tr^ = 2c},jCm — 1 we can write it as a product of two 
matrix product operators each with dimension 2, so that 
is a matrix product operator with dimension at most 
4. This reasoning also holds for models with disorder, 
i.e. where the magnetic field or the couplings depend on 
the lattice site {Bz{j), Jx{j) and Jy{j)), which can not be 
diagonalizcd via Fourier and Bogolubov transformations 



[l5|. Analogous conclusions hold for quasi-free bosonic 
systems. 

This observation demonstrates that a DMRG simula- 
tion in the Heisenberg picture may be considerably more 
efficient, even exact, in cases where the same approach 
in the Schrodinger picture is provably inefficient [lOj. In 
contrast to the Schrodinger picture, the block entangle- 
ment in the Hcisenbcrg picture (considering the four op- 
erators Po, Pi, P2, P3 as basis- vectors of a 4-dim Hilbcrt 
space for each site) is bounded for all times. This dif- 
ference in entanglement scaling in the two pictures obvi- 
ously can not hold for all settings [3. Nonetheless we 
find indications for a saturation in the scaling of block 
entanglement in numerical simulations for more general 
models. 

Numerical results - We now turn to compare the nu- 
merical efficiency of H-DMRG with that of DMRG in 
the Schrodinger picture Q, S S, II ■ For the dynamics of 
pure states, we have seen that there are examples where 
DMRG becomes exact thanks to a very favorable behav- 
ior of entanglement. In general one expects the use of 
H-DMRG to be advantageous only where the entangle- 
ment scaling for the state is drastically worse than for 
the operator to be evolved. This is due to the following 
reason: If a quantum state has the matrix product repre- 

scntation 1*) = ■ ■ ■ 4^] Ih) ® ■ ■ ■ « M 

with matrix dimension d, then the operator has 
the matrix product representation 



...,ijv=0 



(1) 



B 



where P^'^ = A^^^ 



\ii){ii\(E)- . .'S)\iN){iN\, 
(3) 



and hence the B matrices 



have dimension d^ . The matrix product representation of 
an operator is thus expected to require matrix dimension 
d^ in situations, where the representation of a state only 
requires d and is therefore much more efficient. 

The situation is different if decoherence and dissipa- 
tion is present as then the evolution of an operator must 
be considered in both, Heisenberg and Schrodinger pic- 
ture. Dissipation may be described by local Lindblad 
terms leading to a master equation for the dynamics of 
the density matrix g of the form [17| . 



N 



g = -i[H,g]+^^ {2a ^ ga^ - ct+ct^- g - ga+a^ ) 



N 



(4) 



where and r„ are the respective damping rates. When 
the description is transferred into the Heisenberg picture, 




FIG. 1: The time evolution, {a^){t) for a model described by 
eqs. (HJ respectively ^ with parameters A'^ = 10, Bz = 0.8, 
= 0.5, Jy = 0.4, Jz = 0.01, r„ = 0.1 and Ta = 0.1. a: 
The exact solution, b: (5dmrg for DMRG-simulations in the 
Schrodinger picture for rf = 20 (blue), d = 40 (green), d = 60 
(red), d = 80 (cyan) and d — 100 (magenta), c: iSh-dmrg for 
DMRG-simulations in the Heisenberg picture. 



the same dynamics is described by the equation. 



N 



(5) 



for a Heisenberg picture operator X{t). In the following 
we compare the results of numerical simulations in the 
Schrodinger (eq. S]) [11] and Heisenberg picture (eq. O. 

In our first example we choose the parameters of the 
model to be, N = 10, = 0.8, = 0.5, Jy = 0.4, 
Jz = 0.01, Tu = 0.1 and Vd = 0.1 to allow for compari- 
son with exact results. We simulate the time evolution 
of the operator crf(t), where the initial state is all 
spins pointing down in z-direction, |0o) = | 
(^^1 i) = ^1 i))- Figured^ shows the exact solution, 
that is a 2nd order Runge-Kutta integration with time 
steps dt = 0.01 of eq. (g]). AU our DMRG simulations 
also use 2nd order integrations with dt = 0.01. The 
errors of DMRG-simulations in the Schrodinger pic- 
ture, 5dmRG = logio (|((7f)exact(0 - (c^5)DMRG(i)|), 

and the Heisenberg picture, (^h-dmrg = 
logio (l(o-f)exact(O - ('^5)H-DMRG(i)|), ^re shown in 
figures [TJa and c respectively (Calculations with time 
steps dt = 0.001 produced indistinguishable results). 
The Heisenberg picture simulations show a significantly 
higher accuracy than the Schrodinger picture results. 
Moreover the accuracy improvement with increasing 
matrix dimension d is more pronounced in the Heisen- 
berg picture, suggesting an unfavorable scaling of 
entanglement in the Schrodinger picture. 

In a second example we consider a chain of iV = 100 
spins, with slightly different parameters to show the gen- 
erality of our findings, and compare DMRG results in 
the Schrodinger and Heisenberg picture. Here, = 0.8, 
= 0.5, Jy = 0.4, Jz = 0.01, r„ = 0.01 and = 0.1. 




FIG. 2: Dynamics of a chain with N = 100 spins and 
Bz = 0.8, Jo; = 0.5, Jy = 0.4, Jz = 0.01, r„ = 0.01 
and Vd = 0.1. a: {alo)(t) as given by eq. (gj (DMRG) 
for d = 20 (blue), d = 40 (green), d = 60 (red), d = 80 
(cyan) and d = 100 (magenta), b: {(J5o){t) as given by 
eq. © (H-DMRG) for d = 20,40,60,80 and 100. c: 
Wl,o)it)U-[{ala){t)]^2\ as given by eq. O (DMRG) 
for (dl,d2) = (40,20) (blue), (dl,d2) = (60,40) (green), 
(dl,d2) = (80,60) (red) and (dl,d2) = (100,80) (cyan), d: 
- [(«^50>(i)]d2l as given by eq. © (H-DMRG) for 
(dl,d2) = (40, 20), (60, 40), (80, 60) and (100,80). With in- 
creasing bond dimension, results converge much faster in the 
Heisenberg than in the Schrodinger picture. 



Figures [2h and [^Ja show (cr|Q)(t) as calculated in the 
Schrodinger and Heisenberg picture (03) • Since it is 
not possible to compare these values to exact results for 
= 100, we test the convergence of the obtained results 
with increasing matrix dimension, d. This convergence is 
shown in figurel^b for the Schrodinger and in figure[21i for 
the Heisenberg picture, where we plotted the differences 
between results that were obtained with different matrix 

dimensions dl and d2, |[(crfo)(t)]rfi - [(crfo)(t)]rf2l- The 
convergence is found to be much faster in the Heisenberg 
than in the Schrodinger picture. We found the advantage 
of the Heisenberg picture to be of a comparable signif- 
icance as in figure [2] already for a chain with the same 
parameters but = 40 spins. 

In our simulations we use the time evolved block deci- 
mation (TEBD) algorithm introduced in 0]. This algo- 
rithm, at each step, truncates the reduced density ma- 
trices of all considered bipartitions by only keeping the 
states corresponding to their d largest eigenvalues (Here 
d is the dimension of the employed matrices.). We com- 
pute the truncation error by summing up all truncated 
eigenvalues of the reduced density matrices Q at each 
time step. The truncation errors at each time step, et, 
are then cumulatively summed up. The resulting quan- 
tity, e = £(, is thus an upper bound to approximation 
errors due to matrix truncation. 

To enable a comparison of the accuracies of the matrix 
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FIG. 3: Truncation errors e at f = 10 for Schrodinger picture 
DMRG (blue) and H-DMRG (green) for N = 20, 40, 60 and 
80 and fixed bond dimension d = 60. Bz = 0.8, Jx ~ —0.5, 
Jy = 0.4, Jz = 0.1, r„ = 0.1 and Fa = 0.2. H-DMRG 
truncations saturate at e < 0.96. 



truncations, we set Hcisenberg and Schrodinger represen- 
tations on the same footing by normalizing the represen- 
tations of eqs. ^ and ([1]) in the Frobenius norm, i.e. 
for a TO X m matrix X wc set J2Tj=i l-^ijl"^ ~ ^- With 
this normalization, the matrix representations of g and X 
both have the structure of Matrix Product States where 
the eigenvalues of reduced matrices sum up to unity and 
can be interpreted as Schmidt coefRcientsjl3|. 

Since we compare truncation errors in two different 
representations it is not obvious that lower truncation 
errors in one representation imply a better approxima- 
tion for the expectation value of an observable or vice 
versa. Indeed, for the short chains used in the exam- 
ple in figure [1] we found comparable truncation errors in 
both approaches even though the error in the relevant 
observable is much smaller when using H-DMRG [20]. 

On the other hand, the truncation errors appear to 
be significantly lower in H-DMRG for longer chains. We 
have therefore investigated the scaling of the truncation 
error with the system size for fixed bond dimension d. 
Importantly, we observe that in contrast to Schrodinger 
picture DMRG, the truncation errors e are found to 
saturate in H-DMRG for fixed bond dimension if the 
system size is increased. Figure [3] shows the trunca- 
tion errors e at t = 10 for Schrodinger picture DMRG 
(blue) and H-DMRG (green) for d = 60 and chain length 
N = 20,40,60,80 and 100. For the remaining param- 
eters wc chose a set, that is again different from those 
in figures [U and [2] with = 0.8, = -0.5, Jy = 0.4, 

= 0.1, r„ = 0.1 and Pd = 0.2. 

This clear difference in the scaling further corroborates 
the idea that the Heisenberg and Schrodinger picture arc 
qualitatively different in regards to their entanglement 
scaling even beyond the exactly solvable models discussed 
earlier. In fact, this feature hints at a saturation of the 
entanglement of bipartitions in the Heisenberg picture. It 
will be an interesting challenge for future work to provide 
analytical arguments to support the numerical findings 
presented here and to demonstrate more rigorously the 
superior efficiency of H-DMRG when applied to mixed 
state evolutions beyond the numerical findings here. 
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